set.seed(54554)

library(foreign)
library(lpSolve)

D <- read.dta("brader_trim_nomiss.dta")
Y <- D$imm_bin1
M <- ifelse(D$emo > median(D$emo), 1, 0)
X <- D$tone_eth

P000 <- 10000000 * sum(Y==0 & M==0 & X==0)/sum(X==0)
P001 <- 10000000 * sum(Y==0 & M==0 & X==1)/sum(X==1)
P010 <- 10000000 * sum(Y==0 & M==1 & X==0)/sum(X==0)
P011 <- 10000000 * sum(Y==0 & M==1 & X==1)/sum(X==1)
P100 <- 10000000 * sum(Y==1 & M==0 & X==0)/sum(X==0)
P101 <- 10000000 * sum(Y==1 & M==0 & X==1)/sum(X==1)
P110 <- 10000000 * sum(Y==1 & M==1 & X==0)/sum(X==0)
P111 <- 10000000 * sum(Y==1 & M==1 & X==1)/sum(X==1)


f.con <- matrix(c(rep(c(1,1,0,0), 2^3), rep(0, 2^5), #P000
        rep(c(rep(c(1,0),2^2), rep(0,2^3)), 2^2), #P001
        rep(c(rep(c(0,0,1,1), 2^2), rep(0, 2^4)), 2), #P010
        rep(c(0,1,0,1,0,0,0,0), 2^3), #P011
        rep(0, 2^5), rep(c(1,1,0,0), 2^3), #P100
        rep(c(rep(0,2^3), rep(c(1,0), 2^2)), 2^2), #P101
        rep(c(rep(0, 2^4), rep(c(0,0,1,1), 2^2)), 2), #P010
        rep(c(0,0,0,0,0,1,0,1), 2^3), #P011
        rep(1, 2^6) # sum=1
        ), nrow=9, byrow=TRUE)
f.dir <- rep("=", 9)
f.rhs <- c(P000, P001, P010, P011, P100, P101, P110, P111, 10000000)
index <- sample(1:64, 64)
q <- rep(NA, 64)
names(q) <- c("q0000_00","q0000_01","q0000_10","q0000_11",
    "q0001_00","q0001_01","q0001_10","q0001_11",
    "q0010_00","q0010_01","q0010_10","q0010_11",
    "q0011_00","q0011_01","q0011_10","q0011_11",
    "q0100_00","q0100_01","q0100_10","q0100_11",
    "q0101_00","q0101_01","q0101_10","q0101_11",
    "q0110_00","q0110_01","q0110_10","q0110_11",
    "q0111_00","q0111_01","q0111_10","q0111_11",
    "q1000_00","q1000_01","q1000_10","q1000_11",
    "q1001_00","q1001_01","q1001_10","q1001_11",
    "q1010_00","q1010_01","q1010_10","q1010_11",
    "q1011_00","q1011_01","q1011_10","q1011_11",
    "q1100_00","q1100_01","q1100_10","q1100_11",
    "q1101_00","q1101_01","q1101_10","q1101_11",
    "q1110_00","q1110_01","q1110_10","q1110_11",
    "q1111_00","q1111_01","q1111_10","q1111_11")
counter <- 1

for(i in index){
    f.obj <- rep(0, 64)
    f.obj[i] <- 1
    res.max <- lp("max", f.obj, f.con, f.dir, f.rhs)
    res.min <- lp("min", f.obj, f.con, f.dir, f.rhs)
    if (res.max$status == 2 | res.min$status == 2){
        stop("no feasible solution ", counter, "\n")
    }
    if (res.max$objval - res.min$objval < .Machine$double.eps){
        q[i] <- 0
    } else {
        q[i] <- (res.max$objval - res.min$objval)*rbeta(1,1,1) + res.min$objval
    }
    newcons <- rep(0, 64)
    newcons[i] <- 1
    f.con <- rbind(f.con, newcons)
    f.dir <- c(f.dir, "=")
    f.rhs <- c(f.rhs, q[i])
    counter <- counter + 1
}

q <- q/10000000

p <- 0.9 # strength of encouragement
N <- (68+198) * 3000 # Sample Size: original * scalar
probT <- 68/(68+198)

indexfun <- function(j, data){
    index <- as.logical(c(rep(c(rep(0, 2^(j-1)), rep(1, 2^(j-1))), 2^(6-j)),
                          rep(0,ncol(data)-64)))
    return(index)
}

## Generate frequency table for principal strata
Freq  <- round(q*N, 0)
Freq.temp <- round(q*N, 1)

if (sum(Freq)<N) {
    # Note: Rounding errors from multiplication of PS probabilities and sample sizes could produce observation 
    # lengths longer or shorter than what is in the actual sample
    diff <-  N - sum(Freq) #Calculate differences in the sum of frequencies and the sample
    repeat{
        if(diff > 64){
        Freq <- Freq + 1
        diff <- N - sum(Freq)
        } else break
    }
    p.temp <-  Freq.temp - Freq
    thresh  <-  .5 #differences at .4 likely to be under classified
    repeat{
        pot  <-  p.temp >= thresh
        if(sum(pot) < diff) thresh  <-  thresh - .1 #however lower thresholds may be needed
        else break
        cat("Too much rounding down; new threshold at ", thresh, "\n")
    }
    a  <-  sample(names(q)[pot], diff)
    # Update the count of the Frequency of the created data for if clause on top of this loop
    Freq[a]  <-  Freq[a] + 1
}

if (sum(Freq)>N) {
    # Note: Rounding errors from multiplication of PS probabilities and sample sizes could produce observation 
    # lengths longer or shorter than what is in the actual sample
    diff <-  sum(Freq) - N #Calculate differences in the sum of frequencies and the sample
    repeat{
        if(diff > 64){
        Freq <- Freq - 1
        diff <- sum(Freq) - N
        } else break
    }
    p.temp <-  Freq - Freq.temp 
    thresh  <-  .5 #differences at .4 likely to be under classified
    repeat{
        pot  <-  p.temp >= thresh
        if(sum(pot) < diff) thresh  <-  thresh - .1 #however lower thresholds may be needed
        else break
        cat("Too much rounding up; new threshold at ", thresh, "\n")
    }
    a  <-  sample(names(q)[pot], diff)
    # Update the count of the Frequency of the created data for if clause on top of this loop; bc sum(Freq)>N case need to remove cases
    Freq[a]  <-  Freq[a] - 1
}

## Set up fake data frame
data <- as.data.frame(matrix(NA,nrow=N,ncol=64))
colnames(data) <- names(q)

## Fill in PS indicators
# First rep over the rows the number of observations that are in the first stratum
data[,1] <- c(rep(1,Freq[1]),rep(0,N-Freq[1]))
# For the remaining columns this must be done while starting at a different row. 
# e.g. if the first PS had one observation, then the second PS (column) would start at row 2
for(i in 2:64) { 
    a <- sum(Freq[1:i-1])
    data[,i] <- c(rep(0,a),rep(1,Freq[i]),rep(0,N-Freq[i]-a))
}

## Assign T w/ complete randomization
data$T <- sample(c(rep(1,probT*N),rep(0,N-probT*N)), N)

## Assign D w/ block randomization wrt T
data$D <- NA
data$D[data$T==1] <- sample(c(rep(1, probT*N/2), rep(0,probT*N/2)), probT*N)
data$D[data$T==0] <- sample(c(rep(1, (N-probT*N)/2), rep(0, (N-probT*N)/2)), N-probT*N)

## Determine potential M & Y based on PS
index <- indexfun(1,data)
data$M1 <- data$M10 <- apply(data[,index],1,sum)
index <- indexfun(2,data)
data$M0 <- data$M00 <- apply(data[,index],1,sum)
index <- indexfun(3,data)
data$Y11  <-  apply(data[,index],1,sum)
index <- indexfun(4,data)
data$Y10  <-  apply(data[,index],1,sum)  
index <- indexfun(5,data)
data$Y01  <-  apply(data[,index],1,sum)
index <- indexfun(6,data)
data$Y00  <-  apply(data[,index],1,sum)

## Generate potential M for the positively encouraged (Z=1)
Freq.clps <- rep(NA, 16)
names(Freq.clps) <- c("q0000","q0001","q0010","q0011","q0100","q0101","q0110","q0111",
    "q1000","q1001","q1010","q1011","q1100","q1101","q1110","q1111")
for(i in 1:16) Freq.clps[i] <- sum(Freq[(4*(i-1)+1):(4*i)])
WT <- c(.9,.1,.9,.1,.9,.1,.9,.1,.9,.1,.9,.1,.9,.1,.9,.1) # encourage Y(1,1)=0

data$M01 <- data$M11 <- NA
for(y00 in 0:1){
    for(y01 in 0:1){
        for(y10 in 0:1){
            for(y11 in 0:1){
            index <- y00*2^3 + y01*2^2 + y01*2 + y11 + 1
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==1,])>0){
            data$M01[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==1] <- 1
            data$M11[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==1] <- 1
            }
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==0,])>0){
            data$M01[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0] <- 1
            nenc10 <- round(WT[index]*nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0,]),0)
            nunenc10 <- nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0,]) - nenc10
            data$M11[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0] <- sample(c(rep(1,nenc10), rep(0,nunenc10)), nenc10+nunenc10)
            }
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==1,])>0){
            nenc01 <- round(WT[index]*nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1,]),0)
            nunenc01 <- nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1,]) - nenc01
            data$M01[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1] <- sample(c(rep(1,nenc01), rep(0,nunenc01)), nenc01+nunenc01)
            data$M11[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==1] <- 1
            }
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==0,])>0){
            nenc00 <- round(WT[index]*nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==0,]),0)
            nunenc00 <- nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==0,]) - nenc00
            data$M01[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==0] <- sample(c(rep(1,nenc00), rep(0,nunenc00)), nenc00+nunenc00)
            data$M11[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==0] <- data$M01[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==0]
            }
            }
        }
    }
}


## Generate potential M for the negatively encouraged (Z=2)
WT2 <- c(.1,.9,.1,.9,.1,.9,.1,.9,.1,.9,.1,.9,.1,.9,.1,.9) # encourage Y(1,1)=1

data$M02 <- data$M12 <- NA
for(y00 in 0:1){
    for(y01 in 0:1){
        for(y10 in 0:1){
            for(y11 in 0:1){
            index <- y00*2^3 + y01*2^2 + y01*2 + y11 + 1
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==0,])>0){
            data$M02[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==0] <- 0
            data$M12[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==0] <- 0
            }
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==0 & data$M10==1,])>0){
            data$M02[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1] <- 0
            nenc01 <- round(WT2[index]*nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1,]),0)
            nunenc01 <- nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1,]) - nenc01
            data$M12[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==0 & data$M10==1] <- sample(c(rep(0,nenc01), rep(1,nunenc01)), nenc01+nunenc01)
            }
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==0,])>0){
            nenc10 <- round(WT2[index]*nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0,]),0)
            nunenc10 <- nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0,]) - nenc10
            data$M02[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==0] <- sample(c(rep(0,nenc10), rep(1,nunenc10)), nenc10+nunenc10)
            data$M12[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==0] <- 0
            }
            
            if (nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & data$M00==1 & data$M10==1,])>0){
            nenc11 <- round(WT2[index]*nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==1,]),0)
            nunenc11 <- nrow(data[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==1,]) - nenc11
            data$M02[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==1] <- sample(c(rep(0,nenc11), rep(1,nunenc11)), nenc11+nunenc11)
            data$M12[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==1] <- data$M02[data$Y00==y00 & data$Y01==y01 & data$Y10==y10 & data$Y11==y11 & 
                data$M00==1 & data$M10==1]
            }
            }
        }
    }
}

## Assign Z w/ block randomization wrt T
data$Z <- NA
data$Z[data$T==1] <- sample(c(rep(1, probT*N/3), rep(0,probT*N/3), rep(-1,probT*N/3)), probT*N)
data$Z[data$T==0] <- sample(c(rep(1, (N-probT*N)/2), rep(0, (N-probT*N)/2), rep(-1, (N-probT*N)/2)), N-probT*N)

## Determine observed M under parallel for D==0
data$M <- NA
data$M[data$T==1 & data$D==0] <- data$M1[data$T==1 & data$D==0]
data$M[data$T==0 & data$D==0] <- data$M0[data$T==0 & data$D==0]

## Assign M under parallel w/ block randomization wrt T for D==1
nT1D1 <- sum(data$T==1 & data$D==1)
data$M[data$T==1 & data$D==1] <- sample(c(rep(1, nT1D1/2), rep(0, nT1D1/2)), nT1D1)
nT0D1 <- sum(data$T==0 & data$D==1)
data$M[data$T==0 & data$D==1] <- sample(c(rep(1, nT0D1/2), rep(0, nT0D1/2)), nT0D1)

## Determine observed M under encouragement
data$Mz <- NA
data$Mz[data$T==1 & data$Z==0] <- data$M10[data$T==1 & data$Z==0]
data$Mz[data$T==0 & data$Z==0] <- data$M00[data$T==0 & data$Z==0]
data$Mz[data$T==1 & data$Z==1] <- data$M11[data$T==1 & data$Z==1]
data$Mz[data$T==0 & data$Z==1] <- data$M01[data$T==0 & data$Z==1]
data$Mz[data$T==1 & data$Z==-1] <- data$M12[data$T==1 & data$Z==-1]
data$Mz[data$T==0 & data$Z==-1] <- data$M02[data$T==0 & data$Z==-1]

## Determine observed Y under parallel
data$Y <- NA
data$Y[data$T==1 & data$M==1] <- data$Y11[data$T==1 & data$M==1]
data$Y[data$T==1 & data$M==0] <- data$Y10[data$T==1 & data$M==0]
data$Y[data$T==0 & data$M==1] <- data$Y01[data$T==0 & data$M==1]
data$Y[data$T==0 & data$M==0] <- data$Y00[data$T==0 & data$M==0]

## Determine observed Y under encouragement
data$Yz <- NA
data$Yz[data$T==1 & data$Mz==1] <- data$Y11[data$T==1 & data$Mz==1]
data$Yz[data$T==1 & data$Mz==0] <- data$Y10[data$T==1 & data$Mz==0]
data$Yz[data$T==0 & data$Mz==1] <- data$Y01[data$T==0 & data$Mz==1]
data$Yz[data$T==0 & data$Mz==0] <- data$Y00[data$T==0 & data$Mz==0]

## Calculate CME
data$YT1M0 <- data$M0*data$Y11 + (1-data$M0)*data$Y10
data$YT1M1 <- data$M1*data$Y11 + (1-data$M1)*data$Y10
data$YT0M0 <- data$M0*data$Y01 + (1-data$M0)*data$Y00
data$YT0M1 <- data$M1*data$Y01 + (1-data$M1)*data$Y00
data$delta1 <- data$YT1M1 - data$YT1M0
data$delta0 <- data$YT0M1 - data$YT0M0

#Calculate Quantities necessary for computation of bounds
Z111 <- sum(data$Y==1 & data$T==1 & data$M==1 & data$D==1)/sum(data$T==1 & data$M==1 & data$D==1)
Z011 <- sum(data$Y==0 & data$T==1 & data$M==1 & data$D==1)/sum(data$T==1 & data$M==1 & data$D==1)
Z001 <- sum(data$Y==0 & data$T==1 & data$M==0 & data$D==1)/sum(data$T==1 & data$M==0 & data$D==1)
Z110 <- sum(data$Y==1 & data$T==0 & data$M==1 & data$D==1)/sum(data$T==0 & data$M==1 & data$D==1)
Z010 <- sum(data$Y==0 & data$T==0 & data$M==1 & data$D==1)/sum(data$T==0 & data$M==1 & data$D==1)
Z000 <- sum(data$Y==0 & data$T==0 & data$M==0 & data$D==1)/sum(data$T==0 & data$M==0 & data$D==1)
Z101 <- sum(data$Y==1 & data$T==1 & data$M==0 & data$D==1)/sum(data$T==1 & data$M==0 & data$D==1)
Z100 <- sum(data$Y==1 & data$T==0 & data$M==0 & data$D==1)/sum(data$T==0 & data$M==0 & data$D==1)

P111 <- sum(data$Y==1 & data$M==1 & data$T==1 & data$D==0)/sum(data$T==1 & data$D==0)
P101 <- sum(data$Y==1 & data$M==0 & data$T==1 & data$D==0)/sum(data$T==1 & data$D==0)
P001 <- sum(data$Y==0 & data$M==0 & data$T==1 & data$D==0)/sum(data$T==1 & data$D==0)
P011 <- sum(data$Y==0 & data$M==1 & data$T==1 & data$D==0)/sum(data$T==1 & data$D==0)
P110 <- sum(data$Y==1 & data$M==1 & data$T==0 & data$D==0)/sum(data$T==0 & data$D==0)
P100 <- sum(data$Y==1 & data$M==0 & data$T==0 & data$D==0)/sum(data$T==0 & data$D==0)
P000 <- sum(data$Y==0 & data$M==0 & data$T==0 & data$D==0)/sum(data$T==0 & data$D==0)
P010 <- sum(data$Y==0 & data$M==1 & data$T==0 & data$D==0)/sum(data$T==0 & data$D==0)

#Calculate Quantities necessary for computation of bounds
dP000 <- sum(data$Yz==0 & data$Mz==0 & data$T==0 & data$Z==-1)/sum(data$T==0 & data$Z==-1)
dP001 <- sum(data$Yz==0 & data$Mz==0 & data$T==1 & data$Z==-1)/sum(data$T==1 & data$Z==-1)
dP010 <- sum(data$Yz==0 & data$Mz==1 & data$T==0 & data$Z==-1)/sum(data$T==0 & data$Z==-1)
dP011 <- sum(data$Yz==0 & data$Mz==1 & data$T==1 & data$Z==-1)/sum(data$T==1 & data$Z==-1)
dP100 <- sum(data$Yz==1 & data$Mz==0 & data$T==0 & data$Z==-1)/sum(data$T==0 & data$Z==-1)
dP101 <- sum(data$Yz==1 & data$Mz==0 & data$T==1 & data$Z==-1)/sum(data$T==1 & data$Z==-1)
dP110 <- sum(data$Yz==1 & data$Mz==1 & data$T==0 & data$Z==-1)/sum(data$T==0 & data$Z==-1)
dP111 <- sum(data$Yz==1 & data$Mz==1 & data$T==1 & data$Z==-1)/sum(data$T==1 & data$Z==-1)

tP000 <- sum(data$Yz==0 & data$Mz==0 & data$T==0 & data$Z==0)/sum(data$T==0 & data$Z==0)
tP001 <- sum(data$Yz==0 & data$Mz==0 & data$T==1 & data$Z==0)/sum(data$T==1 & data$Z==0)
tP010 <- sum(data$Yz==0 & data$Mz==1 & data$T==0 & data$Z==0)/sum(data$T==0 & data$Z==0)
tP011 <- sum(data$Yz==0 & data$Mz==1 & data$T==1 & data$Z==0)/sum(data$T==1 & data$Z==0)
tP100 <- sum(data$Yz==1 & data$Mz==0 & data$T==0 & data$Z==0)/sum(data$T==0 & data$Z==0)
tP101 <- sum(data$Yz==1 & data$Mz==0 & data$T==1 & data$Z==0)/sum(data$T==1 & data$Z==0)
tP110 <- sum(data$Yz==1 & data$Mz==1 & data$T==0 & data$Z==0)/sum(data$T==0 & data$Z==0)
tP111 <- sum(data$Yz==1 & data$Mz==1 & data$T==1 & data$Z==0)/sum(data$T==1 & data$Z==0)

sP000 <- sum(data$Yz==0 & data$Mz==0 & data$T==0 & data$Z==1)/sum(data$T==0 & data$Z==1)
sP001 <- sum(data$Yz==0 & data$Mz==0 & data$T==1 & data$Z==1)/sum(data$T==1 & data$Z==1)
sP010 <- sum(data$Yz==0 & data$Mz==1 & data$T==0 & data$Z==1)/sum(data$T==0 & data$Z==1)
sP011 <- sum(data$Yz==0 & data$Mz==1 & data$T==1 & data$Z==1)/sum(data$T==1 & data$Z==1)
sP100 <- sum(data$Yz==1 & data$Mz==0 & data$T==0 & data$Z==1)/sum(data$T==0 & data$Z==1)
sP101 <- sum(data$Yz==1 & data$Mz==0 & data$T==1 & data$Z==1)/sum(data$T==1 & data$Z==1)
sP110 <- sum(data$Yz==1 & data$Mz==1 & data$T==0 & data$Z==1)/sum(data$T==0 & data$Z==1)
sP111 <- sum(data$Yz==1 & data$Mz==1 & data$T==1 & data$Z==1)/sum(data$T==1 & data$Z==1)

P0 <- c(dP000,dP010,dP100,dP110,tP000,tP010,tP100,tP110,sP000,sP010,sP100,sP110)
P1 <- c(dP001,dP011,dP101,dP111,tP001,tP011,tP101,tP111,sP001,sP011,sP101,sP111)
Q0 <- c(dP010+dP110,tP010+tP110,sP010+sP110)
Q1 <- c(dP011+dP111,tP011+tP111,sP011+sP111)

#Single Experiment
l.delta.t <- max(-P001-P011, -P000-P001-P100, -P011-P010-P110)
u.delta.t <- min(P101+P111, P000+P100+P101, P010+P110+P111)
l.delta.c <- max(-P100-P110, -P001-P101-P100, -P011-P111-P110)
u.delta.c <- min(P000+P010, P011+P111+P010, P000+P001+P101)
S.t <- cbind(l.delta.t,u.delta.t)
S.c <- cbind(l.delta.c,u.delta.c)

#Parallel Design
l.delta.t <- max(-P001-P011, -P011-P010-P110-P001+Z001, -P000-P001-P100-P011+Z011, -P001-P011+Z001-Z111, -P001+P101-Z101, -P011+P111-Z111)
u.delta.t <- min(P101+P111, P010+P110+P101+P111-Z101, P000+P100+P101+P111-Z111, P101+P111+Z001-Z111, P111-P011+Z011, P101-P001+Z001)
l.delta.c <- max(-P100-P110, -P011-P111-P110-P100+Z000, -P001-P101-P100-P110+Z110,-P100-P110+Z100-Z010, -P110+P010-Z010, -P100+P000-Z100)
u.delta.c <- min(P000+P010, P011+P111+P010+P000-Z100, P000+P001+P101+P010-Z010, P000+P010+Z100-Z010, P010-P110+Z110, P000-P100+Z100)
P.t <- cbind(l.delta.t,u.delta.t)
P.c <- cbind(l.delta.c,u.delta.c)

#Encouragement Design
bounds.bidir.cmpl <- function(P,Q,dir){
    f.obj <- c(rep(0,16), 0,0,1,0, 0,0,1,0, 0,-1,0,0, 0,-1,0,0, 
            0,0,-1,0, 0,0,-1,0, 0,1,0,0, 0,1,0,0, rep(0,16))
    f.con <- matrix(c(
            rep(c(1,1,1,0),8), rep(0,32), #dP00t
            rep(c(rep(c(0,0,0,1),4), rep(0,16)),2), #dP01t
            rep(0,32), rep(c(1,1,1,0),8), #dP10t
            rep(c(rep(0,16), rep(c(0,0,0,1),4)),2), #dP11t
            rep(c(1,1,0,0),8), rep(0,32), #tP00t
            rep(c(rep(c(0,0,1,1),4), rep(0,16)),2), #tP01t
            rep(0,32), rep(c(1,1,0,0),8), #tP10t
            rep(c(rep(0,16), rep(c(0,0,1,1),4)),2), #tP11t
            rep(c(1,0,0,0),8), rep(0,32), #sP00t
            rep(c(rep(c(0,1,1,1),4), rep(0,16)),2), #sP01t
            rep(0,32), rep(c(1,0,0,0),8), #sP10t
            rep(c(rep(0,16), rep(c(0,1,1,1),4)),2), #sP11t
            rep(c(rep(0,12), rep(1,4)),4), #dQ
            rep(c(rep(0,8), rep(1,8)),4), #tQ
            rep(c(rep(0,4), rep(1,12)),4), #sQ
            rep(1,64) #sum=1
            ), nrow=16, byrow=T)
    f.dir <- rep("=", 16)
    f.rhs <- c(P,Q,1)
    r <- lp(dir, f.obj, f.con, f.dir, f.rhs)
    r$objval
}

bounds.bidir.popl <- function(P,Q,dir){
    f.obj <- c(rep(0,16), 0,0,1,1, 0,0,1,1, -1,-1,0,0, -1,-1,0,0, 
            0,0,-1,-1, 0,0,-1,-1, 1,1,0,0, 1,1,0,0, rep(0,16))
    f.con <- matrix(c(
            rep(c(1,1,1,0),8), rep(0,32), #dP00t
            rep(c(rep(c(0,0,0,1),4), rep(0,16)),2), #dP01t
            rep(0,32), rep(c(1,1,1,0),8), #dP10t
            rep(c(rep(0,16), rep(c(0,0,0,1),4)),2), #dP11t
            rep(c(1,1,0,0),8), rep(0,32), #tP00t
            rep(c(rep(c(0,0,1,1),4), rep(0,16)),2), #tP01t
            rep(0,32), rep(c(1,1,0,0),8), #tP10t
            rep(c(rep(0,16), rep(c(0,0,1,1),4)),2), #tP11t
            rep(c(1,0,0,0),8), rep(0,32), #sP00t
            rep(c(rep(c(0,1,1,1),4), rep(0,16)),2), #sP01t
            rep(0,32), rep(c(1,0,0,0),8), #sP10t
            rep(c(rep(0,16), rep(c(0,1,1,1),4)),2), #sP11t
            rep(c(rep(0,12), rep(1,4)),4), #dQ
            rep(c(rep(0,8), rep(1,8)),4), #tQ
            rep(c(rep(0,4), rep(1,12)),4), #sQ
            rep(1,64) #sum=1
            ), nrow=16, byrow=T)
    f.dir <- rep("=", 16)
    f.rhs <- c(P,Q,1)
    r <- lp(dir, f.obj, f.con, f.dir, f.rhs)
    r$objval
}

BE.t <- c(bounds.bidir.popl(P1, Q0, "min"), bounds.bidir.popl(P1, Q0, "max"))
BE.c <- c(-bounds.bidir.popl(P0 ,Q1, "max"), -bounds.bidir.popl(P0, Q1, "min"))
num.BET.t.lo <- bounds.bidir.cmpl(P1, Q0, "min")
num.BET.t.up <- bounds.bidir.cmpl(P1, Q0, "max")
num.BET.c.lo <- -bounds.bidir.cmpl(P0, Q1, "max")
num.BET.c.up <- -bounds.bidir.cmpl(P0, Q1, "min")
denom.BET.t <- P1[1] + P1[3] - P1[9] - P1[11]
denom.BET.c <- P0[1] + P0[3] - P0[9] - P0[11]
BET.t <- c(num.BET.t.lo/denom.BET.t, num.BET.t.up/denom.BET.t)
BET.c <- c(num.BET.c.lo/denom.BET.c, num.BET.c.up/denom.BET.c)


# actual proportion encouraged
p <- mean(data$M11 != data$M10 | data$M01 != data$M00)

delta1 <- mean(data$delta1)
delta0 <- mean(data$delta0)
tdelta1 <- mean(data$delta1[(data$M12==0 & data$M10==0 & data$M11==1) | 
                                        (data$M12==0 & data$M10==1 & data$M11==1)])
tdelta0 <- mean(data$delta0[(data$M02==0 & data$M00==0 & data$M01==1) | 
                                        (data$M02==0 & data$M00==1 & data$M01==1)])


# Figure 2
pdf("fig2.pdf", width=12, height=4)
par(mfrow=c(1,2), mar=c(5.1,1.1,1.1,1.1), oma=c(0,0,3,0), cex=0.8)

plot(0,0,type="n",xlim=c(-1,1),ylim=c(-.5,3.5), main="", ylab="", cex.lab=1.5,
    xlab="", axes = FALSE, bty="n")
axis(1, at = seq(-1,1,by=.2))
axis(3, at = c(-1,0,1), labels = FALSE)
points(rep(delta1,3), c(3,2,1), pch=16, cex=1.5)
points(tdelta1, 0, pch=16, cex=1.5)
lines(S.t, c(3,3), lwd=3)
lines(P.t, c(2,2), lwd=3)
lines(BE.t, c(1,1), lwd=3)
lines(BET.t, c(0,0), lwd=3)
points(S.t, c(3,3), pch="|")
points(P.t, c(2,2), pch="|")
points(BE.t, c(1,1), pch="|")
points(BET.t, c(0,0), pch="|")
abline(v=0)
abline(v=-1)
abline(v=-.8, lty="dotted")
abline(v=-.6, lty="dotted")
abline(v=-.4, lty="dotted")
abline(v=-.2, lty="dotted")
abline(v=1)
abline(v=.8, lty="dotted")
abline(v=.6, lty="dotted")
abline(v=.4, lty="dotted")
abline(v=.2, lty="dotted")
text(-0.95, 3, "Single Experiment", pos=4)
text(-0.95, 2, "Parallel         ", pos=4)
text(-0.95, .5, "Encouragement     ", pos=4)
text(delta1, 3.25, expression(paste("Population Effect ", bar(delta), "(",1,")")))
text(delta1, 2.25, expression(paste("                            ", bar(delta), "(",1,")")))
text(delta1, 1.25, expression(paste("                            ", bar(delta), "(",1,")")))
text(tdelta1, .25, expression(paste("Complier Effect ", bar(delta)^" *", "(",1,")")))

plot(0,0,type="n",xlim=c(-1,1),ylim=c(-.5,3.5), main="", ylab="", cex.lab=1.5,
    xlab="", axes = FALSE, bty="n")
axis(1, at = seq(-1,1,by=.2))
axis(3, at = c(-1,0,1), labels = FALSE)
points(rep(delta0,3), c(3,2,1), pch=16, cex=1.5)
points(tdelta0, 0, pch=16, cex=1.5)
lines(S.c, c(3,3), lwd=3)
lines(P.c, c(2,2), lwd=3)
lines(BE.c, c(1,1), lwd=3)
lines(BET.c, c(0,0), lwd=3)
points(S.c, c(3,3), pch="|")
points(P.c, c(2,2), pch="|")
points(BE.c, c(1,1), pch="|")
points(BET.c, c(0,0), pch="|")
abline(v=0)
abline(v=-1)
abline(v=-.8, lty="dotted")
abline(v=-.6, lty="dotted")
abline(v=-.4, lty="dotted")
abline(v=-.2, lty="dotted")
abline(v=1)
abline(v=.8, lty="dotted")
abline(v=.6, lty="dotted")
abline(v=.4, lty="dotted")
abline(v=.2, lty="dotted")
text(0.95, 3, "Single Experiment", pos=2)
text(0.95, 2, "         Parallel", pos=2)
text(0.95, .5, "Encouragement", pos=2)
text(delta0, 3.25, expression(paste("Population Effect ", bar(delta), "(",0,")")))
text(delta0, 2.25, expression(paste("                            ", bar(delta), "(",0,")")))
text(delta0, 1.25, expression(paste("                            ", bar(delta), "(",0,")")))
text(tdelta0, .25, expression(paste("Complier Effect ", bar(delta)^" *", "(",0,")")))

title(main="Average Indirect Effects", cex.main=1.5, outer=T)

dev.off()
